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Prediction  of  coastal  processes,  including  waves,  currents,  and  sediment  transport,  can  be  obtained  from  a 
variety  of  detailed  geophysical-process  models  with  many  simulations  showing  significant  skill.  This 
capability  supports  a  wide  range  of  research  and  applied  efforts  that  can  benefit  from  accurate  numerical 
predictions.  However,  the  predictions  are  only  as  accurate  as  the  data  used  to  drive  the  models  and,  given  the 
large  temporal  and  spatial  variability  of  the  surf  zone,  inaccuracies  in  data  are  unavoidable  such  that  useful 
predictions  require  corresponding  estimates  of  uncertainty.  We  demonstrate  how  a  Bayesian-network  model 
can  be  used  to  provide  accurate  predictions  of  wave-height  evolution  in  the  surf  zone  given  very  sparse  and/or 
inaccurate  boundary-condition  data.  The  approach  is  based  on  a  formal  treatment  of  a  data-assimilation 
problem  that  takes  advantage  of  significant  reduction  of  the  dimensionality  of  the  model  system.  We 
demonstrate  that  predictions  of  a  detailed  geophysical  model  of  the  wave  evolution  are  reproduced  accurately 
using  a  Bayesian  approach.  In  this  surf-zone  application,  forward  prediction  skill  was  83%,  and  uncertainties  in 
the  model  inputs  were  accurately  transferred  to  uncertainty  in  output  variables.  We  also  demonstrate  that  if 
modeling  uncertainties  were  not  conveyed  to  the  Bayesian  network  (i.e.,  perfect  data  or  model  were 
assumed),  then  overly  optimistic  prediction  uncertainties  were  computed.  More  consistent  predictions  and 
uncertainties  were  obtained  by  including  model-parameter  errors  as  a  source  of  input  uncertainty.  Improved 
predictions  (skill  of  90%)  were  achieved  because  the  Bayesian  network  simultaneously  estimated  optimal 
parameters  while  predicting  wave  heights. 

Published  by  Elsevier  B.V. 


1.  Introduction 

The  coastal  environment  is  characterized  by  extreme  variability. 
Physical  variables  that  exhibit  significant  spatial  and  temporal 
evolution  include  waves,  currents,  and  bathymetry.  This  is  particularly 
true  in  the  surf  zone,  where  waves  transition  from  non-breaking  to 
breaking  conditions,  transferring  momentum  to  drive  currents, 
sediment  transport,  and  bathymetric  change.  Numerous  numerical 
geophysical  models  are  available  that  make  predictions  of  these 
processes.  These  process  models  can  also  be  inherently  statistical, 
wherein  only  quantities  that  are  formally  averaged  over  several  wave 
periods  (and,  often,  averaged  over  the  water  column)  are  simulated. 
For  instance,  temporally  averaged  statistical  properties  of  waves  can 
be  accurately  predicted  by  SWAN  (Simulating  Waves  Nearshore,  Booij 
et  al„  1999;  Ris  et  al.,  1999)  given  accurate  bathymetry,  water  levels, 
description  of  the  frequency-directional  spectra  at  boundary  condi¬ 
tions,  and  specification  of  a  number  of  tuning  parameters.  Similarly, 
wave-averaged  currents  can  be  predicted  by  models  such  as  Delft-3D 
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(Lesser  et  al.,  2004;  Reniers  et  al.,  2007)  or  ADC1RC  (Westerink  et  al., 
2008)  that  have  clearly  demonstrated  predictive  skill.  At  higher 
resolution,  wave-resolving  Boussinesq  models  also  have  excellent 
predictive  skill  (Chen  et  al.,  2003)  and  an  ability  to  simulate  at  the 
shorter  time  scales  given  correspondingly  high-resolution  time-series 
data  on  the  model  boundaries.  This  increased  fidelity  may  be  necessary 
for  predicting  sediment  transport  and  bathymetric  evolution  (Hen¬ 
derson  et  al.,  2004).  However,  the  higher  temporal  resolution  of 
Boussinesq  models  comes  with  increased  computational  cost  along 
with  the  more  demanding  specification  of  boundary  conditions. 

At  the  present  level  of  hydrodynamic  modeling  capability, 
improved  predictive  skill  typically  depends  on  improving  the 
accuracy  of  model-boundary  conditions  (Plant  et  al.,  2009),  rather 
than  on  refinements  of  the  model  parameterizations,  indicating  that 
the  geophysical  theory  governing  nearshore  processes  is  relatively 
mature.  Modeling  of  sediment-transport  processes  is  an  exception  to 
this  statement,  where  skillful  predictions  can  depend  strongly  on 
choice  of  model  parameterization  (Ruessink  et  al„  2007).  Parameter 
dependence  is  primarily  due  to  the  use  of  wave-averaged  models  that 
do  not  resolve  physically  important  processes  such  as  the  details  of 
the  wave  boundary  layer  and  higher-than-second-order  moments  of 
the  near-bed  velocities  (Henderson  et  al.,  2004;  Hsu  et  al.,  2006). 
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These  neglected  details  must  be  absorbed  into  the  available  free 
parameters.  The  bottom  line  is  that  the  predictive  skill  of  the  present 
modeling  capability  will  largely  depend  on  uncertainties  in  model 
inputs,  model  parameters,  or  both. 

Finally,  real-world  implementation  of  coastal-process  models, 
useful  to  emergency  managers,  the  military,  lifeguards,  and  beach- 
goers  who  need  to  make  decisions  based  on  rapid  environmental 
forecasts,  is  becoming  a  necessity.  In  this  scenario,  the  fidelity  of 
boundary  and  initial  conditions  required  by  numerical  models  may  be 
so  poor  that  the  inadequacies  of  the  input  data  can  only  be  improved 
by  additional  (and  often  costly)  observations  of  processes  in  the 
interior  of  the  model  domain.  This  sets  up  a  complex  situation  where 
optimal  solutions  require  the  use  of  some  sort  of  assimilation  model  to 

(1)  map  incomplete  and  inaccurate  data  to  the  boundary  conditions, 

(2)  correct  model  errors  in  the  interior  of  the  domain,  and  (3)  extend 
intrinsic  model  outputs  (e.g.,  wave  height)  to  related  observables 
(e.g.,  video  or  radar  observations,  Bell,  1999;  van  Dongeren  et  al., 
2008)  that  can  contribute  to  model  improvement.  Presently,  formal 
data-assimilation  schemes  that  have  been  applied  to  nearshore 
coastal  processes  have  focused  primarily  on  scientific  evaluation  of 
model  physics  (Feddersen  et  al.,  2004),  studying  the  sensitivity  of 
inverse  solutions  to  the  scale  of  variability  (Kurapov  et  al.,  2007),  and 
evaluating  unknown  model  parameters  (Plant  et  al.,  2004;  Ruessink  et 
al.,  2007),  rather  than  to  solving  the  operational  problem  of 
inadequate  boundary  conditions. 

The  lack  of  regularly  applied  data  assimilation  to  nearshore  coastal 
models  is  likely  due  to  computational  constraints  of  the  state-of-the- 
art  models,  which  often  perform  simulations  at  computational  times 
that  exceed  the  simulation  time.  Furthermore,  the  coupled  wave, 
current,  and  sediment-transport  problem  is  nonlinear  and  requires 
significant  model  iteration  to  converge  toward  the  best  inverse 
solutions  (Kurapov  et  al.,  2007).  To  get  an  idea  for  the  dimensionality 
of  the  problem,  consider  a  coastal  region  that  might  span  100  km2, 
with  a  model  resolution  of  1 00  m2,  which  requires  evaluation  at 
1  million  grid  positions.  At  each  position,  there  may  be  roughly  10 
variables  of  interest  (e.g.,  depth,  wave  height,  energy  dissipation,  etc.). 
The  models  might  have  a  time  step  of  1  s,  and  a  forecast  would 
estimate  hourly  averaged  conditions  requiring  evaluation  of  a  total  of 
1010  values  to  describe  the  coastal  environment.  A  typical  assimilation 
scheme  might  require  10  iterations,  pushing  the  total  bookkeeping 
load  to  1011  values  per  coastal-process  forecast. 

Although  formal  data  assimilation  is  possible,  its  implementation 
is  difficult  in  coastal  and  surf-zone  environments  because  (1)  it  is 
computationally  intensive;  (2)  it  requires  specification  of  often 
unknown  quantities  (such  as  the  model  error);  (3)  it  is  unwieldy  in 
the  face  of  a  large  number  of  uncertainties;  and  (4)  it  is  unmanageable 
in  the  face  of  a  large  number  of  observations.  Nonetheless,  in  practical 
applications  of  coastal  prediction,  forecasts  will  be  expected  even 
from  boundary  and  forcing  data  with  substandard  accuracy  (e.g., 
outdated  bathymetric  survey  and  hydrodynamic  observations,  or 
forecasts  with  significant  uncertainty).  At  a  minimum,  model 
initialization  will  require  some  form  of  data  assimilation  to  correct 
the  imperfect  model  inputs.  For  instance,  data  interpolation  of  one 
form  or  another  is  a  commonly  used  (and  abused)  form  of  data 
assimilation  (Ooyama,  1987;  Plant  et  al„  2002).  Regardless  of 
approach,  practical  applications  can  be  supported  only  if  the  data- 
assimilation  method  offers  some  reduction  in  observation  and  model 
errors  and,  most  importantly,  provides  a  useful  estimate  of  prediction 
uncertainties. 

In  this  two-part  paper,  we  present  a  new  application  of  Bayesian 
statistical  modeling  to  surf-zone  modeling  and  assimilation.  Part  I  is 
devoted  to  describing  the  modeling  framework  and  applying  it  to 
several  case  studies  focused  on  forward  modeling.  Part  II  is  devoted  to 
inverse  problems  and  extending  the  framework  to  a  veiy  flexible  set 
of  assimilation  problems.  We  continue  with  Part  I  in  Section  2  (Model 
formulation)  by  reviewing  Bayes  Rule  and  then  developing  a  specific 


application  for  surf-zone  wave  modeling.  In  Section  3  (Application), 
we  describe  a  data  set  used  for  training  the  Bayesian  network  and 
provide  hindcast  and  forecast  prediction  examples.  These  examples 
highlight  the  role  that  data  and  model  uncertainty  play  in  affecting 
prediction  errors.  Lastly,  in  Section  4  (Discussion),  we  explore  the 
problem  of  obtaining  required  uncertainty  estimates  and  demonstrate 
their  impact  on  model  skill. 

2.  Model  formulation 

2.1.  Review  of  Bayes  Rule 

In  contrast  to  variational  data  assimilation,  mentioned  in  the 
previous  section,  we  suggest  that  a  Bayesian-inference  approach  can 
be  used  to  combine  available  measurements  with  nearshore-process 
models  to  make  statistically  robust  forecasts.  The  Bayesian  approach  is 
formally  consistent  with  sensitivity-based  variational  assimilation 
(Wilde  and  Berliner,  2007).  Bayes  Rule  is 

p(Fi|Oj)  =p(0J\Fi)p(Fi)/p(0j),  (1) 

where  the  left-hand  side  of  Eq.  (1)  is  the  updated  conditional 
probability  of  a  particular  forecast,  Fit  given  a  particular  set  of 
observations,  Oj.  The  forecast  might  include  both  initial  and  boundary 
conditions.  The  observations  might  be  obtained  near  the  boundaries 
or  in  the  interior  of  the  model  domain.  Variational  data  assimilation 
strives  to  return  the  conditional  mean  value  of  the  solution:  Fj  =  S,  p 
(F,  |  Oj)  Fj.  A  Bayesian  approach  is  more  general  and  strives  to  estimate 
the  likelihood  of  many  possible  forecasts.  The  first  term  on  the  right- 
hand  side  of  Eq.  (1)  is  the  inverse  of  the  left-hand  side  and  is  the 
likelihood  of  the  observations  if  the  forecast  is  known.  This  term  can 
include  both  model  and  observation  errors.  That  is,  if  the  model  and 
measurements  were  error  free,  then  an  observation  would  be  likely 
only  if  it  equaled  the  forecast  value.  In  reality,  there  are  numerous 
errors  causing  spread  in  the  likelihood  function. 

The  next  term  on  the  right-hand  side  of  Eq.  (1)  is  the  prior 
probability  of  each  forecast.  This  is  what  is  known  about  the  problem 
before  new  data  are  available.  It  might  be  the  result  of  a  prior 
assimilation  cycle  or  derived  from  climatology.  The  mean  value 
obtained  from  this  distribution  is  equivalent  to  the  “background"  or 
“best-guess”  solution  used  in  variational  ocean  data  assimilation 
(Benett,  2002)  or  in  optimal  interpolation  (Ooyama,  1987).  Finally, 
the  last  term  is  a  normalization  factor  to  account  for  the  total 
likelihood  of  the  observations.  In  variational  data  assimilation  and 
optimal  interpolation,  this  term  is  solved  by  inverting  the  data 
covariance  and  is  often  responsible  for  large  computational  costs. 
Here,  it  is  estimated  through  integration  over  all  forecast  possibilities. 

A  primary  advantage  to  the  Bayesian  approach  to  data  assimilation 
is  that  the  probability  distribution  of  the  forecast  is  estimated.  This 
allows  for  both  data  and  forecast  to  have  non-Gaussian  probability 
distributions,  which  may  be  crucial  to  describing  strongly  nonlinear 
nearshore  processes.  In  principle,  this  means  that  the  Bayesian 
approach  will  lead  to  more  accurate  forecast  statistics,  including 
estimates  of  the  most  likely  forecast  and  its  uncertainty.  A  possible 
disadvantage  of  the  Bayesian  approach  is  that  it  may  increase  the 
dimensionality  of  the  problem.  If  each  F;  is  a  single  forecast  with  1010 
quantities  to  describe  it  (i.e.,  the  typical  numerical  model  system 
mentioned  earlier),  then  the  complete  Bayesian  approach  requires 
that  we  track  the  joint  probability  of  these  variables  against 
observations,  which  might  also  be  any  of  the  1010  quantities.  The 
problem  would  have  exploded  to  a  dimension  of  10lo°! 

in  practice,  however,  the  utility  of  the  Bayesian  approach  is  most 
evident  when  the  problem  of  interest  can  be  boiled  down  to  a  much- 
reduced  dimensionality.  (This  statement  actually  applies  to  most 
other  data-assimilation  schemes,  which,  for  instance,  do  not  attempt 
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to  compute  the  model-data  sensitivity  at  full  resolution,  Kurapov 
et  al.,  2007.)  Here,  we  will  apply  three  standard  approaches  to 
reducing  the  problem's  dimensionality.  First,  we  forecast  along  a  1- 
dimensional  (1-D)  slice  of  the  spatial  domain,  extending  from  the 
shoreline  across  the  width  of  the  surf  zone  into  offshore  water  depths 
where  waves  do  not  usually  break.  Second,  we  consider  only  a  subset 
of  the  entire  spatial  domain  by  focusing  on  a  relatively  few  locations  of 
interest.  Thirdly,  we  will  consider  only  a  subset  of  the  possible 
processes  in  the  surf-zone  region  that  we  deem  most  important, 
specifically  wave  shoaling  and  breaking. 

2.2.  A  surf-zone  Bayesian  model 

To  test  our  hypothesis  that  a  Bayesian  approach  to  data 
assimilation  can  be  exploited  in  the  surf  zone,  we  selected  the 
wave-evolution  model  described  by  Thornton  and  Guza  (1983) 
(hereinafter  referred  to  as  TG83),  which  assumes  Rayleigh-distribut¬ 
ed  wave  heights  throughout  the  model  domain.  Our  intent  is  not  to 
substantiate  a  particular  approach  over  another— since  numerous 
refinements  to  this  model  have  been  developed— but  to  emphasize 
that  variations  in  implementation  do  not  change  the  essential 
elements  of  this  data-assimilation  example.  The  inherently  1-D 
TG83  model  predicts  the  evolution  of  the  root  mean  square  (rms) 
wave  height  given  estimates  at  the  offshore  boundary  of  the  wave 
height,  wave  direction,  and  wave  period.  The  model  does  not  resolve 
the  different  directions  or  frequencies  that  might  contribute  to  a 
spectrum  of  waves,  because  waves  are  assumed  to  be  narrowly 
distributed  around  a  peak  period  and  direction.  Thus,  we  can  write  the 
forecast  from  this  model  at  any  location,  xk,  as, 

Fi(.xk)  =  {HJ,.  =fiinct.({h0,hu...,hk,H0,a0J,y,B}i),  (2) 

where  h  is  the  water  depth,  H0  is  the  offshore  rms  wave  height,  and  a 0 
is  the  offshore  peak-wave  direction.  Variables  with  subscript  k  are 
spatially  varying.  Additional  parameters  are  a  peak-wave  period  (T),  a 
critical  wave-breaking  criterion  (y),  and  a  breaking  efficiency  (B).  All 
of  the  inputs  can  be  considered  random  variables,  because  even  the 
model  parameters  must  be  estimated  from  the  data  (Ruessink  et  al., 
2003).  Therefore,  the  ith  solution  is  obtained  from  a  particular  choice 
of  inputs.  Fig.  1  shows  an  example  of  the  spatially  extensive  (in  1 
horizontal  dimension)  input  and  outputs  of  this  model. 

The  waves  shoal  offshore  from  the  outer  sandbar  and  then  begin  to 
break,  decreasing  in  height  until  they  reach  the  shoreline.  The  change 


cross-shore  distance  (m) 

Fig.  1.  Wave  height  (dashed  line)  simulation  over  measured  bathymetry  (solid  line). 
Vertical  lines  mark  reduced  Bayesian  network's  sample  locations. 


in  the  wave-height  profile  will  depend  on  the  offshore  wave  height, 
period,  direction,  and  the  underlying  bathymetry.  Using  the  numerical 
TG83  model,  it  is  not  possible  to  make  a  prediction  unless  all  of  these 
variables  are  known.  In  particular  for  our  case,  the  bathymetry  must 
be  known  at  approximately  500  locations  in  order  to  provide  1-m 
model  resolution  across  the  domain.  The  dimensionality  of  the  full 
model  is,  4*500  +  3,  which  includes  the  spatially  varying  fields  of 
input  bathymetry  and  output  wave  height,  direction,  and  dissipation, 
plus  the  wave  period  and  two  model  parameters. 

2.3.  Reduced  model 

We  expect  that  the  utility  of  wave-model  forecasts  does  not 
necessarily  rely  on  providing  all  of  the  output  predictions  possible 
from  the  model  (O'Reilly  and  Guza,  1998).  For  instance,  accurate 
prediction  of  the  wave  height  very  near  the  beach  might  be  necessary, 
rather  than  at  all  locations  in  the  model  domain.  Or,  if  a  higher- 
resolution  model  is  to  be  initialized  from  a  low-resolution  model,  only 
the  data  at  the  boundary  of  the  high-resolution  model  are  required. 
From  this  point  of  view,  the  unused  model  details  are  irrelevant.  This 
implies  that  our  data-assimilation  problem  should  only  retain  a 
forecast  at  a  small  subset  of  locations.  For  the  surf-zone  example,  we 
choose  to  retain  information  at  our  model's  boundary,  at  an 
intermediate  location,  and  at  one  location  close  to  the  shoreline.  We 
have,  in  an  ad  hoc  way,  just  reduced  the  spatial  dimensionality  from 
500  locations  to  only  3. 

The  Bayesian  versions  of  this  reduced  model  inputs  and  outputs 
are 

and  (3) 

Oj  =  {h0,H0.h,.HvH2.h2J.a0,y,B}j. 

The  spatial  locations  are  given  by  subscripts  (0  for  the  offshore 
boundary,  1  for  the  intermediate  position,  and  2  nearshore).  A 
variable  with  a  circumflex  indicates  that  it  is  an  observation  (which 
can  be  derived  from  any  number  of  sources,  including  other  models) 
and  the  others  are  predictions.  We  have  included  the  model 
parameters  y  and  B  as  both  observations  and  predictions.  However, 
without  calibration  data,  these  values  will  be  obtained  from  a  prior 
distribution.  Initially,  we  choose  to  constrain  these  parameters  to 
y  =  0.34  and  B  =  0.8  (Haines  and  Sallenger,  1994)  in  order  to  focus 
attention  on  the  remaining  variables.  We  will  loosen  this  constraint 
when  we  discuss  parameter  sensitivity  later. 

A  schematic  diagram  of  the  Bayesian  network  representing  the 
reduced  model  is  shown  in  Fig.  2.  This  follows  a  hierarchical  approach 
(Ihlera  et  al.,  2007),  which  considers  only  a  subset  of  the  total  joint 
correlations  that  are  possible  for  this  system.  The  network  variables 
are  grouped  both  spatially  (from  left  to  right)  and  by  process  (top  to 
bottom).  The  arrows  connecting  specific  variables  or  groups  of 
variables  represent  joint  correlations  between  those  variables.  These 
correlations  will  be  estimated  from  measurements  or  from  model 
predictions  as  part  of  the  training  or  calibration  process.  Although  it  is 
entirely  possible  to  represent  correlations  between  each  variable  and 
all  other  variables,  we  have  chosen  to  take  advantage  of  a  further 
reduction  in  model  dimensionality  by  representing  only  a  few  of  the 
most  important  correlations  (each  is  labeled  with  a  number).  The 
choice  of  which  correlations  to  retain  is  guided  by  our  knowledge  of 
the  wave-evolution  processes.  For  example,  we  know  that  the  wave 
height  at  location  X!  depends  on  the  wave  height  at  x0  (correlation 
#5),  and  additional  relationships  are  expected  for  the  water  depths  at 
x0  and  Xj  (correlations  #7  and  #8),  wave  direction  and  period 
(correlation  #2),  and  the  model  parameters  (correlation  #4). 

In  addition,  we  know  that  other  processes  contribute  to  correla¬ 
tions  that  would  not  be  represented  by  the  wave  model,  but  would  be 
expressed  in  the  observations.  For  instance,  wave  height,  direction, 
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Fig.  2.  Surf-zone  model  representation  as  a  Bayesian  network. 


and  period  can  be  correlated  to  each  other  (#1  and  #2)  due  to 
constraints  associated  with  fetch  and  storm  characteristics.  Water 
depths  are  correlated  to  each  other  (#11  and  #12)  due  to  mass 
conservation  and  coherent  variations  associated  with  sandbars. 

The  Bayesian  network  of  the  reduced  model  in  the  assimilation 
problem  has  several  useful  properties,  including  that  it 

(1 )  it  achieves  a  necessary  reduction  in  dimensionality; 

(2)  it  focuses  a  detailed  forecast  model  on  variables  and  locations 
of  interest; 

(3)  it  estimates  forecast  uncertainty; 

(4)  it  can  be  used  to  inversely  estimate  boundary  conditions  for  a 
detailed  forecast  model  (to  be  shown  in  Part  II);  and 

(5)  it  can  assimilate  diverse  inputs,  including  depth,  wave  height, 
wave  direction,  and  dissipation  (and  more,  as  we  will  also  show 
in  Part  II). 

Nonetheless,  the  Bayesian  network  still  relies  on  detailed 
numerical  estimates  and/or  observations  because  the  joint  probability 
that  relates  p(F  |  0)  to  p(0  |  F),  p(F),  and  p(0)  must  be  estimated 
through  training  based  on,  essentially,  Monte  Carlo  simulation.  The 
proposed  efficiency  in  the  Bayesian  approach  is  in  the  fact  that  once  all 
relevant  scenarios  have  been  encountered,  the  detailed  numerical 
model  used  in  the  “training"  process  can  be  discarded.  This  is 
particularly  powerful  in  the  nearshore  environment  where  dynamics 
periodically  return  to  similar  beach  states  (Plant  et  al„  2006;  Wright  et 
al.,  1985)  and  even  to  identical  bathymetric  profiles  at  regular 
intervals  (Ruessink  and  Kroon,  1994;  Wijnberg  and  Terwindt,  1995). 
Rather  than  re-computing  the  same  inputs  with  a  detailed  model  to 
get  nearly  the  same  outputs,  the  Bayesian  approach  stores  not  only 
the  functional  mapping  between  input  and  output,  but  also  the 
uncertainty  in  the  output  due  to  uncertainty  in  the  input. 

3.  Application 

3.1.  Bayesian-network  discretization 

We  constructed  a  representation  of  Eq.  (1)  and  Eq.  (2)  with  the 
reduced  model  Eq.  (3)  using  a  numerical  implementation  of  a  Bayesian 
network  called  Netica  (Norsys,  1990-2007).  We  have  already  justified 
the  spatial  simplifications  of  the  reduced  model.  An  additional  level  of 
simplification  is  required  to  represent  the  continuous  physical 


processes  with  a  small  number  of  discrete  states.  That  is,  a  continuous 
range  of  natural  wave  heights,  H,  will  be  represented  with  a  finite 
distribution  of  values,  Hn,  ranging  from  some  minimum  value  (e.g., 
0  m)  to  a  reasonable  maximum  value  (e.g.,  the  maximum  observed 
value  or  a  fraction  of  deepest  water  depth,  since  we  know  that  wave 
height  will  be  limited  to  this  by  breaking).  It  is  not  necessary  to  have 
constant  intervals  (Ffn  +  1  —  Hn^Hn  —  Ffn„i).  The  design  of  this 
discretization  requires  balancing  a  number  of  competing  interests. 

(1)  The  intervals  should  be  as  wide  as  possible  to  minimize 
computational  effort. 

(2)  The  intervals  should  be  narrow  enough  to  resolve  the  expected 
uncertainty  and  provide  forecast  utility.  For  instance,  if  the 
expected  rms  prediction  error  were  crH,  then  the  interval 
should  be  no  larger  than  aH. 

(3)  The  intervals  need  to  be  wide  enough  to  collect  at  least  a  few 
“hits”  from  the  available  data  and  Monte  Carlo  simulations.  The 
objective  is  to  approximate  the  joint  probabilities  of  interest  by, 
essentially,  compiling  histograms.  If  an  interval  contains  only  a 
small  number  of  samples,  the  histogram  is  not  well  con¬ 
strained.  In  practice,  we  seed  all  of  the  joint  probabilities  with  a 
uniform  distribution  so  that  poorly  constrained  intervals  will, 
at  worst,  return  a  uniform  distribution  as  a  result. 

Constraint  #2  suggests  that  the  network  resolution  should  be  as 
detailed  as  possible  while  supporting  constraints  #1  and  #3.  This  is 
not  necessarily  the  case.  There  are  intrinsic  errors,  such  as  model 
error,  that  can  be  easily  captured  by  preventing  the  intervals  from 
shrinking  unnecessarily.  Thus,  constraint  #2  can  be  interpreted  to 
mean  that  the  intervals  should  be  no  finer  than  about  aH  as  well. 
Additionally,  this  constraint  can  be  interpreted  from  a  practical  point 
of  view,  and  the  interval  should  be  no  finer  than  what  is  necessary  for 
the  appropriate  interpretation  of  the  results.  For  example,  an 
application  to  beach  safety  might  require  only  estimating  whether 
the  wave  height  exceeds  a  critical  value.  Then,  only  two  bins  are 
required,  split  at  the  critical  value.  We  adhere  to  these  design 
guidelines  using  an  adaptive  method,  described  in  the  next  section. 

3.2.  Training  data  set 

The  proposed  Bayesian  network  could  be  developed  for  any 
arbitrary  location.  We  selected  the  Army  Corps  of  Engineers  Field 
Research  Facility  (FRF),  in  Duck,  NC.  At  this  site  there  is  a  long  time 
series  of  bathymetry  and  wave  observations.  These  data  are  sufficient 
to  drive  a  forward  model  of  wave  shoaling  and  breaking 
corresponding  to  a  large  range  of  prediction  scenarios.  We  emphasize 
that  locations  with  less  extensive  data  sets  could  also  work  within  this 
modeling  concept,  but  that  the  source  data  from  this  location  are  well 
established  and  easy  to  obtain.  Fig.  3  shows  histograms  of  daily- 
averaged  data  for  the  period  1980  to  1996,  using  essentially  the  same 
data  set  as  that  analyzed  by  Plant  et  al.  (1999).  The  wave  data  were 
measured  approximately  every  3  h.  The  bathymetry  was  measured 
approximately  monthly  and  was  interpolated  to  a  daily  interval  using 
an  objective  interpolation  scheme.  The  computational  domain  of  the 
TG83  mode  extended  offshore  to  the  8-m  contour  where  wave  data 
were  measured.  The  measured  bathymetry  was  spatially  interpolated 
to  1-m  intervals  across  the  domain.  The  measured  inputs  and  model 
parameters  (y  and  B)  were  inputted  to  the  TG83  model,  and  the 
output  wave-height  predictions  were  stored  at  offshore,  intermediate 
and  onshore  locations.  Note  that  the  purpose  of  the  modeling  exercise 
at  this  stage  is  to  obtain  data  to  be  used  to  estimate  correlations  both 
among  and  between  model  inputs  and  outputs.  Although  this  could  be 
accomplished  with  synthetic  data,  the  field  data  set  provides  accurate 
prior  probabilities  for  the  input  and  output  variables.  The  interpola¬ 
tion  scheme  produced  smooth  results  in  both  space  and  time  and 
filtered  out  some  short  spatial  and  temporal  variability.  A  potential 
impact  of  this  is  that  some  wave  conditions  may  be  paired  with  some 


N.G.  Plant,  K.T.  Holland  /  Coastal  Engineering  58  (20 11 )  1 19-130 


123 


profile  conditions  that  are  not  fully  consistent  with  nature.  If  this 
problem  is  significant,  new  observations  may  be  identified  as  being 
inconsistent  (highly  improbable)  or,  more  likely,  this  simply  adds  to 
predicted  uncertainty. 

Using  the  observations  and  corresponding  simulations,  we  initially 
assigned  minimum  bin  widths  to  each  variable.  These  widths  were 
crh  =  0.25  m,  cth  =  0.2  m,  oT=2  s,aa  =  15°.  The  bins  were  then  altered 
such  that  each  would  contain  20%  of  the  data  (e.g.,  five  bins  if  the  data 
were  uniformly  distributed).  If  this  initial  scheme  produced  bins  that 
were  narrower  than  the  minimum  value,  the  bin  width  was  increased 
such  that  there  might  be  fewer  than  five  bins  for  a  variable.  The 
histograms  shown  in  Fig.  3  reflect  the  resolution  that  resulted  using 
this  scheme. 

The  observation/simulation  results  were  read  into  Netica,  which 
assimilated  these  data  into  tables  representing  the  joint  probability  of 
the  wave  model  and  input  data.  Initially,  the  joint  probabilities  were 
given  a  uniform  distribution.  Then,  these  distributions  were  updated 
to  yield  the  maximum  likelihood  probabilities,  given  the  available 
data.  That  is,  the  data  were  used  to  solve  an  inverse  model  for  the 
Bayesian  network  itself,  treating  each  joint  probability  as  an  unknown 
parameter.  There  were  15,435  unknown  joint  probabilities.  These 
were  constrained  with  39,037  observed  and  simulated  values.  Since 
we  assume  that  the  probabilities  of  adjacent  intervals  are  somewhat 
correlated,  there  is  no  requirement  that  the  number  of  inputs  exceeds 
the  number  of  unknown  network  variables.  And,  we  expect  that  many 
scenarios  are  unlikely  (for  instance,  shallow  depths  at  all  locations 
would  not  occur),  such  that  the  initial  uniform  distribution  adequately 
initializes  these  unobserved  states.  Furthermore,  it  is  possible  to 
assimilate  data  with  missing  values  for  some  of  the  variables.  For 
instance,  there  were  numerous  occasions  where  the  wave  direction  at 
the  offshore  boundary  was  unknown.  These  cases  simply  distribute 
their  probability  over  all  possible  values  of  wave  direction  while 


constraining  the  joint  probability  of  the  remaining  non-missing 
variables. 

3.3.  Prediction  example 

Once  the  Bayesian  network  is  “trained"  on  the  available  observa¬ 
tions  and  simulations,  we  can  use  it  to  make  forward  and  inverse 
predictions  (the  latter  will  be  described  in  Part  II).  The  purpose  of  the 
forward  prediction  is  to  demonstrate  the  accuracy  of  the  reduced 
model  relative  to  the  original  model  as  well  as  to  illustrate  how 
degraded  or  inaccurate  input  data  affect  the  uncertainty  of  the  output 
prediction. 

3.3.1.  Prior  prediction 

The  first  example  of  the  Bayesian-network  prediction  has  already 
been  demonstrated  in  Fig.  3.  It  is  the  prior  prediction,  which  is  available 
before  we  constrain  and  update  the  network  with  specific  case  data.  In 
this  situation,  the  prior  prediction  captures  our  understanding  of  the 
likelihood  of  particular  outcomes  based  on  climatology.  For  instance, 
some  variables  actually  change  very  little  (e.g.,  the  water  depth  at  the 
seaward  boundary),  whereas  others  are  extremely  variable  (the  wave 
height  at  the  seaward  boundary).  Data  will  be  compared  to  this  prior 
set  of  distributions  via  Eq.  (1 ).  The  immediate  value  of  the  prior  is  to 
identify  the  consistency  of  the  data.  Data  that  are  consistent  with  the 
prior  will  lead  to  increased  certainty  for  all  of  the  variable's 
distributions.  Data  that  are  inconsistent  will  result  in  very  uncertain 
distributions.  Thus,  the  prior  is  useful  for  evaluating  the  quality  of  the 
data  as  well  as  for  updating  the  predictions. 

3.3.2.  Forward  prediction 

The  next  example  of  the  Bayesian  network  (Fig.  4)  demonstrates  its 
ability  to  reproduce  specific  predictions,  given  data  describing 
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Fig.  4.  Comparison  of  Bayesian  network  prior  prediction  (red  histograms)  and  updated 
prediction  (cross-hatched  histogram)  to  (a)  detailed  input  bathymetry  (blue  lines)  and 
(b)  detailed  wave  shoaling  and  breaking  predictions  (green  lines). 


boundary  conditions.  Using  the  Bayesian  network,  we  reproduce  the 
simulations  shown  in  Fig.  4  with  offshore  wave  height  known  to  be  in 
the  range  of  0.4  to  0.8  m,  period  of  6  to  10  s,  unknown  direction,  and 
water-depth  distributions  centered  on  5.5,  3.8,  and  1.5  m  at  the 
offshore,  intermediate,  and  onshore  locations.  Given  this  information, 
the  predicted  wave  height  (shown  without  bathymetry  in  Fig.  4b)  at 
the  intermediate  location  is  most  likely  to  be  in  the  range  of  0.4  to 
0.6  m  (60%  likelihood),  but  there  is  some  likelihood  for  values  in  the 
range  of  0.6  to  0.8  m  (40%).  At  the  onshore  location,  there  is  more 
confidence  for  the  0.4  to  0.6  m  range  (80%),  with  the  remaining  20% 
likelihood  falling  in  the  0.6  to  0.8  m  range.  In  all  cases,  the  prediction 
that  was  obtained  by  updating  the  Bayesian  network  (black  histo¬ 
grams)  has  a  significantly  narrower  distribution  than  the  prior 
distribution  (red  histograms),  and  the  updated  predictions  are,  as 
expected,  more  consistent  with  the  results  from  the  TG83  wave  model. 

The  next  forward  prediction  (Fig.  5)  replaces  the  well-constrained 
input  wave  height  at  the  offshore  boundary  with  more  uncertain 
information.  The  input  wave  height  was  updated  with  an  input 
probability  distribution  that  was  uniformly  distributed  over  the  range 
of  0  to  1  m.  This  is  meant  to  illustrate  the  effects  of  measurement  error 
or,  alternatively,  forecast  uncertainty  if  one  was  initializing  the  wave 
model  with  output  from  a  numerical  weather  prediction.  This 
example  also  indicates  the  utility  of  having  a  so-called  “informative” 
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Fig.  5.  Comparison  of  Bayesian  prediction  with  increased  offshore  wave  height 
uncertainty. 


prior  distribution.  Even  though  the  offshore  boundary  data  were 
inputted  with  uniform  distributions,  the  updated  probability  dis¬ 
tributions  (Fig.  5)  of  the  offshore  wave  height  are  very  similar  to  the 
priors.  The  Bayesian  network  “knows”  that  moderate  wave  heights 
are  more  likely  than  either  extremely  high  or  low  values  and  corrects 
the  inaccurate  input  data. 

The  result  of  the  increased  input  uncertainty  is  the  increased 
prediction  uncertainty  at  the  landward  locations.  Specifically,  the 
wave-height  prediction  at  the  landward  location  has  degraded  to  just 
50%  likelihood  for  values  in  the  0.4-  to  0.6-m  range,  compared  to  80% 
likelihood  in  the  previous  example.  While  the  prediction  uncertainty 
at  the  landward-most  location  has  increased,  it  is  clear  that  the  0.4-  to 
0.6-m  estimate  is  far  more  likely  than  any  other  value.  This  result  is  due 
to  the  wave-breaking  processes  that  act  to  reduce  the  wave-height 
variability  close  to  shore.  Thus,  uncertainty  in  some  input  variables 
does  not  necessarily  have  the  same  impact  on  all  output  variables. 

The  next  prediction  example  (Fig.  6)  illustrates  the  Bayesian 
networks'  ability  to  identify  inconsistent  data.  In  this  case,  the 
offshore  wave  conditions  used  in  the  first  example  were  applied,  but 
unlikely  depth  observations  were  provided  as  input  (depths  at  all 
three  locations  were  set  to  the  shallowest  levels  with  100%  certainty). 
The  result  is  somewhat  degraded  certainty  of  the  wave  height  at  the 
intermediate  location  (i.e.,  roughly  50%  likelihood  for  both  0.4-  to  0.6- 
and  0.6-  to  0.8-m  levels).  However,  at  the  onshore  location,  the  wave- 
height  prediction  probability  is  distributed  uniformly  over  all 


Fig.  6.  Prediction  with  inconsistent  depth  inputs. 
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possibilities,  which  is  worse  than  the  uncertainty  associated  with  the 
prior  prediction.  This  increase  in  uncertainty  indicates  that  this 
specific  input  scenario  is  one  that  has  no  prior  information.  Either  the 
Bayesian  network  must  be  trained  to  encompass  this  new  (perhaps 
unrealistic)  input,  or  it  demands  that  the  data  be  inspected  for  errors. 

3.3.3.  Field  evaluation 

The  previous  examples  compared  the  Bayesian  network  to  the 
model  output  and  data  that  were  used  to  train  it.  This  hindcast  simply 
demonstrated  how  the  Bayesian  network  responds  to  different  types 
of  input  and  demonstrates  that  it  can  recover  the  detailed  TG83  model 
prediction.  Next,  we  compare  the  Bayesian  predictions  to  new  field 
observations  of  surf-zone  wave  heights.  The  Bayesian  network  was 
not  trained  with  these  new  surf-zone  measurements,  relying  on  the 
previous  training  as  before. 

The  new  field  data  were  obtained  as  part  of  the  Ducl<94  experiment 
(Birkemeier  and  Thornton,  1994).  We  used  the  data  from  17  pressure 
sensors  (Fig.  7)  that  measured  wave  heights  (Elgar  et  al„  1997; 
Gallagher  et  al,  1998).  Bathymetry  data  were  obtained  by  interpola¬ 
tion  of  daily  surveys.  Tides  observed  on  a  nearby  pier  were  used  to 
correct  the  bathymetry  to  obtain  depth  estimates  appropriate  to  each 
observed  wave  height.  As  before,  the  offshore  wave  conditions  and 
depths  at  all  locations  were  used  to  update  the  Bayesian  network  and 
obtain  forward  predictions  of  the  wave  heights  at  the  intermediate  and 
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landward  locations.  We  focus  the  remainder  of  our  analysis  on  the 
Bayesian  prediction  of  the  wave  height  at  the  onshore  location  only. 

We  assume,  initially,  that  all  of  the  data  are  error  free— at  least  to 
the  level  of  precision  required  to  specify  the  Bayesian  network's 
inputs.  Fig.  8  shows  a  comparison  of  the  observed  and  predicted  wave 
heights  at  the  onshore  location.  The  predictions  track  the  observed 
wave  height  variations  that  are  associated  with  changes  in  water 
depth  (tide)  and  offshore  wave  height.  The  skill  of  the  Bayesian 
network's  mean  prediction  is  79%  with  a  mean  error  of  —  0.04  m  (i.e., 
an  under-prediction)  and  rms  error  of  0.15  m.  (The  skill  was  defined 
as  1  —  o£/o%,  where  mse  is  the  mean-square  error  between 
Bayesian  predictions  and  the  observations  and  is  the  variance  of 
the  observed  nearshore  wave  height.)  The  most  likely  prediction  has 
about  the  same  mean  and  rms  errors  and  skill  ( —  0.05  m,  0.14  m,  and 
80%,  respectively).  However,  it  is  clear  from  Fig.  8  that  the  prediction 
skill  is  not  uniformly  accurate.  During  the  second  half  of  the  time 
series  (when  the  bathymetric  surveys  indicated  strong  alongshore 
variability),  the  prediction  systematically  underestimated  the  obser¬ 
vations.  In  fact,  a  more  serious  error  lies  in  the  predicted  uncertainty 
(indicated  by  the  width  of  the  shadings),  which  does  not  include  the 
observations  and,  hence,  is  under-predicted  during  the  period  with 
the  highest  waves. 

Calculation  of  wave  heights  using  the  original  TG83  wave  model  for 
this  interval  shows  that  the  tendency  for  under-prediction  of  the 
highest  wave  heights  is  actually  an  intrinsic  numerical  modeling  error 
and  is  not  due  to  the  Bayesian  re-formulation.  We  will  demonstrate 
that  this  systematic  error  can  be  reduced  by  allowing  for  uncertainty  in 
the  model  parameters,  which  have  been  held  constant  so  far.  Fig.  9 
shows  a  comparison  of  observed  and  predicted  wave  heights  at  the 
onshore  prediction  location.  The  mean  prediction  error  of  the  wave 
model  (labeled TG83)  was  —  0.03  m;  the  rms  errorwas0.07  m,and  the 
prediction  skill  was  97%.  The  improved  skill  and  roughly  7  cm  of 
reduced  rms  error  were  obtained  at  the  expense  of  providing  the  entire 
bathymetric  profile  to  the  TG83  model— rather  than  just  the 
bathymetry  at  the  three  locations  required  by  the  reduced  Bayesian 
network— and  numerically  solving  a  differential  equation  for  each  new 
set  of  observations  (no  matter  how  slightly  they  may  have  changed). 

The  high-resolution  numerical  model  cannot  do  anything  to 
identify  systematic  prediction  errors  (such  as  errors  associated  with 
the  high  waves  at  the  end  of  the  study  period).  On  the  other  hand,  we 
demand  that  the  Bayesian  network  predicts  both  best  estimate  (e.g., 
mean  and  most  likely)  values  and  uncertainties.  The  point  of  tracking 
uncertainty  is  to  alert  users  of  the  predictions  to  situations  when  there 
is  a  high  probability  of  prediction  error.  Unfortunately,  for  the  second 
half  of  the  data  set,  there  is  a  serious  failure  in  the  Bayesian  prediction 
to  identify  prediction  uncertainty  because  the  Bayesian  network  was, 
naively,  based  on  an  inaccurate  wave  model. 

4.  Discussion 

4.1.  Quantifying  Bayesian-prediction  skill 

The  Bayesian-assimilation  approach  turns  the  typical  model- 
evaluation  procedure  on  its  head.  Typically,  a  model  is  calibrated  to 
make  its  best  prediction  through  tuning  of  parameters  in  order  to  fit  a 
hindcast  prediction  through  a  cloud  of  observations.  The  assumption 
is  that  the  model-data  mismatch  will  be  used  to  evaluate  the 
probability  distribution  of  tunable  parameters  and  prediction  errors. 
The  Bayesian-network  prediction  already  includes  the  probability 
distributions.  Thus,  it  is  possible  to  evaluate  the  network  with  each 
independent  observation.  This  is  done  by  evaluating  the  likelihood  of 
an  observation,  given  the  updated  prediction.  The  likelihood  is 


Fig.  7.  Instrument  locations  and  bathymetry  on  (a)  04  and  (b)  24  October  1994. 


F:  =  o\ 


(4) 


126 


N.G.  Plant,  I(.T.  Holland  /  Coastal  Engineering  58  (201 1 )  1 19-130 


a 


date 


b 


date 


Fig.  8.  Observation  of  (a)  offshore  wave  height  and  (b)  comparison  of  observations  and  predictions  at  the  onshore  location.  The  shading  indicates  the  predicted  uncertainty  with 
lightest  shading  for  the  95%  confidence  interval,  medium  shading  for  90%  confidence,  and  dark  for  the  50%  confidence.  The  blue  dots  are  the  observed  surf-zone  values,  and  magenta 
circles  indicate  the  most  likely  prediction  (i.e.,  the  Bayesian-network  bin  that  contains  most  probability  at  each  sample  time).  The  yellow  crosses  mark  the  mean  prediction  (i.e.,  the 
expected  value  from  the  Bayesian  prediction). 


where  0,  is  the  subset  of  observables  that  are  made  available  to  the 
Bayesian  network  for  prediction  and  O'j  is  the  independent  observa¬ 
tion  that  was  withheld  from  the  prediction.  For  instance,  0,  might 
include  the  bathymetry  at  all  locations  as  well  as  the  wave  height  at 
the  offshore  location,  and  Oj  might  be  the  wave  height  at  the  onshore 
location. 


An  approach  to  testing  the  Bayesian-network  prediction  is  to 
compare  it  to  a  competing  model.  One  such  competing  model  is  the 
prior  probability,  which  is  determined  before  the  network  is  updated 
with  specific  observations.  If  the  likelihood  of  the  data  is  increased 
under  the  updated  network,  then  it  will  have  a  positive  log  likelihood 
ratio: 


0.5  1 

Hrms  predicted  (m) 


LRj  =  log  jp  (f;  [Oj)  f  =  Q,  I  —  log  {p(Fi)f|  =  ct }  ■  (5) 

A  likelihood  ratio  exceeding  1  indicates  a  significant  improvement 
in  the  prediction.  This  will  result  if  the  Bayesian  network  confidently 
predicts  a  distribution  centered  on  the  observation.  A  ratio  less  than  0 
is  worse  than  the  prior  prediction.  This  can  result  in  the  case  of  a 
confident  prediction  that  misses  the  verification  data  (such  as  the 
high-wave  prediction  errors).  Thus,  the  likelihood  ratio  simulta¬ 
neously  tests  the  skill  of  the  best  prediction  (which  could  be  the  mean 
or  mode  of  the  distribution)  and  tests  the  skill  of  the  uncertainty 
prediction  by  penalizing  both  overly  certain  and  uselessly  uncertain 
predictions. 

The  likelihood  ratio  summed  over  all  observations  and  predictions 
that  are  shown  in  Fig.  8  was  —  1000.  That  is,  even  though  the  mean 
predictions  were  very  skillful,  the  estimate  of  the  uncertainty  was  not 
accurate.  The  uncertainties  were  overly  confident  such  that  the 
observations  were  more  likely  to  come  from  the  prior  network 
prediction  (which  was  the  same  for  all  sample  times)  than  the 
updated  prediction.  One  explanation  for  this  is  the  omission  of 
information  about  input  errors.  Particularly  important  are  depth 


Fig.  9.  Comparison  of  predicted  and  observed  wave  heights. 
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Fig.  10.  Comparison  of  Duck94  observations  to  predictions  for  input  bathymetric  uncertainty  equal  to  (a)  0,  (b)  1,  and  (c)  4  times  the  true  bathymetric  uncertainty  values.  Including 
wave-parameter  uncertainty  (d)  improves  the  prediction  and  leads  to  more  consistent  estimates  of  the  uncertainty. 
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errors  and  model  errors.  Depth  errors  may  result  from  survey 
inaccuracies  and  alongshore  variability.  Since  the  underlying  process 
model  assumed  alongshore  uniform  bathymetry,  an  obvious  solution 
is  to  include  alongshore  variability  as  an  additional  measure  of 
uncertainty.  Model  errors  might  result  from  using  the  wrong  model- 
parameter  values. 

4.2.  Prediction  sensitivity  to  bathymetric  uncertainty 

To  test  whether  including  input  uncertainty  could  improve  the 
likelihood  score,  we  generated  input  distributions  for  the  depth  data 
from  the  computed  residual  errors  between  survey  data  near  the 
Ducl<94  instrument  array  and  the  interpolated  profile  used  as  input  for 
both  the  Bayesian  and  TG83  models.  Because  the  alongshore  spacing  of 
survey  data  was  between  25  and  50  m,  residuals  were  computed  from 
the  data  within  100  m  of  the  alongshore  position  of  the  array.  This 
selection  of  data  is  somewhat  arbitrary,  but  it  balanced  a  need  for 
sufficient  samples  to  allow  calculation  of  a  meaningful  residual  error, 
yet  still  remains  focused  on  the  location  where  the  waves  were 
measured.  Fig.  10  shows  the  updated  predictions  that  utilize  the  input 
uncertainty,  which  was  assumed  to  be  normally  distributed  with 
standard  deviations  that  were,  on  average,  0.65,  0.15,  0.12  m  at  the 
three  input  locations.  The  uncertainty  in  the  Bayesian  network's 
predictions  changed  most  in  the  second  half  of  the  time  series,  which 
corresponds  to  the  development  of  a  crescentic-bar  shape  (Fig.  7).  The 
likelihood  ratio  increased  to  LR=  300  when  bathymetric  uncertainty 
was  included,  indicating  more  accurate  uncertainty  estimation.  The 
error  statistics  of  the  Bayesian  mean  prediction  were  unchanged. 

Since  the  local  alongshore  variability  might  not  represent  all  of  the 
impacts  of  bathymetry  errors  and  it  certainly  does  not  represent  the 
impact  of  intrinsic  modeling  errors,  it  is  possible  to  use  bathymetric 
uncertainty  as  a  tuning  parameter  to  optimize  the  prediction  accuracy. 
We  progressively  increased  the  amount  of  uncertainty  input  to  the 
Bayesian  network  by  factors  of  two  and  found  that  the  highest 
likelihood  ratio  (i.e.,  the  most  improvement  over  the  prior  prediction) 
was  achieved  for  input  uncertainty  equal  to  4  times  the  actual  values. 
The  best  ratio  was  LR=  1800.  This  result  indicates  that,  perhaps,  the 
missing  cross-shore  detail  in  the  bathymetric  variability  (due  to  using 
a  subset  of  spatial  locations)  plus  the  missing  alongshore  detail  (due  to 
the  assumption  of  alongshore  uniformity)  should  be  included  together 
in  the  input  uncertainty.  Alternatively,  the  source  of  prediction  error 
lies  elsewhere  (e.g.,  with  the  model  parameters)  and  the  use  of 
unrealistically  high  depth  uncertainty  is  simply  a  cover-up  for  the  real 
problem.  In  fact,  it  is  clear  that  a  side  effect  of  improving  the  accuracy  of 
the  Bayesian  estimate  of  uncertainty  is  to  reduce  the  variability  of  the 
predicted  response.  For  instance,  in  Fig.  10c,  the  mean  and  most  likely 
predictions  do  not  exhibit  a  strong  tidal  modulation  for  the  first  half  of 
the  record.  And,  the  Bayesian  network  continues  to  under-predict  the 
wave  height  in  the  second  half  of  the  record. 

Table  1  summarizes  the  error  statistics  of  the  most  landward 
location  for  all  the  updating  methods  described  so  far.  Three 
additional  methods  are  included.  One,  labeled  “Onshore  Hrms”,  used 
the  measured  surf-zone  wave  height  at  the  onshore  prediction 
location  as  input  to  the  Bayesian  network.  That  is,  we  showed  the 
Bayesian  network  the  answer.  The  purpose  of  this  test  is  to  identify 
the  degraded  skill  due  to  the  resolution  of  the  probability  bins.  As  the 
bins  became  coarser,  the  network  would  not  adequately  resolve  the 
best  answer  even  when  that  answer  was  known.  In  fact,  there  was 
very  little  degradation  of  the  model  skill  when  compared  to  the  TG83 
results:  skill  reduction  from  97%  to  94%.  This  test  is  also  a  good 
reference  since  it  indicates  the  maximum  possible  likelihood  score 
(  =  5683). 

The  method  labeled  “No  Bathymetry”  did  not  include  any 
bathymetry  data  as  input  to  the  Bayesian  prediction— that  is,  only 
the  prior  information  from  the  training  on  the  historical  data  set  was 
utilized.  Only  the  offshore  wave  observations  were  used  to  make 


Table  1 

Error  and  skill  statistics  for  mean  wave  height  prediction  at  the  onshore  location. 


Update  method 

Mean  error 

RMS  error 

Skill 

Likelihood 

(m) 

(m) 

(%) 

ratio 

TG83  model 

-0.03 

0.07 

97 

n/a 

No  uncertainty 
Bathymetric 

-0.04 

0.15 

79 

-1011 

uncertainty 
Uncertainty  =  lcr 

-0.04 

0.15 

83 

320 

Uncertainty  =  2cr 

-0.04 

0.15 

83 

1386 

Uncertainty  =  4cr 

-0.04 

0.16 

81 

1859 

Uncertainty  =  8cr 

-0.03 

0.16 

77 

1815 

No  Bathymetry 

7  uncertainty 

-0.06 

0.20 

68 

1483 

Uncertainty  =  0.05 

0.02 

0.12 

89 

2400 

Uncertainty  =  0.1 

0.03 

0.12 

87 

2338 

Uncertainty  =  0.4 

0.08 

0.14 

87 

2003 

Onshore  Hrms 

-0.00 

0.07 

94 

5683 

No  Offshore  Hrms 

-0.13 

0.30 

12 

-157 

Prior 

-0.13 

0.30 

0 

0 

predictions  in  this  case.  The  method  labeled  “No  Offshore  Hrms”  used 
the  bathymetiy  observations  but  did  not  use  the  measured  wave 
parameters.  Mean  and  rms  errors  were  highest  for  the  case  where  the 
wave-height  observations  were  omitted,  indicating  that  observations 
of  bathymetry  were  of  secondary  importance.  With  no  update  in 
bathymetry,  the  network  propagates  the  known  bathymetric  uncer¬ 
tainty  to  the  wave-height  prediction,  scoring  well  (1483)  in  the 
likelihood-ratio  test.  If  the  wave  height,  direction,  and  period  are  not 
updated,  the  “No  Offshore  Hrms”  scenario  does  not  score  well  in  the 
likelihood-ratio  test  (  —  157)  and  has  poor  skill  (12%).  An  outcome  of 
this  comparison  is  that  we  can  place  relative  values  on  updated 
bathymetry  vs.  updated  wave  observations.  In  this  case,  the  offshore 
wave-height  data  contribute  about  six  times  more  than  the 
bathymetry  data  to  skill  improvement.  This  result  depends,  of  course, 
on  the  duration  of  the  field  test,  which  was  long  enough  to  experience 
a  wide  range  of  wave  conditions,  but  shorter  than  needed  to 
experience  the  full  range  of  bathymetric  states. 

4.3.  Prediction  sensitivity  to  wave-model  parameter  uncertainty 

Alternatively,  we  could  address  the  prediction  errors  by  consid¬ 
ering  model  errors.  The  most  straightforward  error  source  is  probably 
the  choice  of  the  two  free  parameters  (y  and  B )  that  control  the  wave¬ 
breaking  process  (Fig.  2).  Others  have  shown  that  the  optimal  values 
of  these  parameters  actually  depend  on  the  incident  wave  conditions 
and  water  depth  (Apotsos  et  al.,  2008;  Ruessink  et  al„  2003).  To 
account  for  model-parameter  uncertainty,  we  generated  new  simula¬ 
tions  by  varying  the  wave-breaking  parameters  uniformly  (y  =  0.17 
to  0.68  and  B  =  0.4  to  1.6  in  increments  that  were  10%  of  the  original 
values),  choosing  one  combination  of  parameters  at  a  time  and  re¬ 
running  the  wave  model  for  each  of  the  points  in  the  20-year  time 
series.  This  resulted  in  a  total  of  1,334,017  simulations  that  were  then 
used  to  re-train  a  Bayesian  network,  modified  to  include  the  model 
parameters  as  inputs. 

Then,  the  Ducl<94  data  set  was  used  to  update  and  evaluate  the 
Bayesian-network  predictions.  Because  it  has  been  shown  that  the 
optimal  values  of  the  two  model  parameters  are  correlated  (Roelvink 
and  Broker,  1993),  we  chose  to  only  allow  uncertainty  in  the  breaking 
parameter,  y,  and  we  constrained  the  other  parameter  at  the  original 
value  (B  =  0.8).  The  Bayesian  network  has,  presumably,  learned  the 
correlation  between  the  model  parameters;  it  should  be  able  to  cope 
with  this  constraint.  The  resulting  skill,  error,  and  likelihoods  from 
this  sensitivity  study  are  shown  in  Table  1.  The  best  prediction  of  the 
nearshore  wave  height  was  obtained  when  the  uncertainty  in  y  was 
equal  to  0.05  (1  standard  deviation),  or  about  10%  of  the  original 
parameter  value  (0.34).  The  mean  and  rms  wave-height  errors  were 
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Fig.  11.  Probability  distribution  (contours)  and  mean  value  (heavy  solid  line)  of 
Bayesian  estimate  of  the  wave-breaking  parameter  value  (y)  vs.  offshore  wave  height. 

lower,  and  the  skill  and  likelihood  ratios  were  higher  than  any  other 
choice  of  uncertainty  values  applied  to  y.  We  tried  some  scenarios 
where  both  bathymetric  and  y  uncertainties  were  allowed  (not 
shown),  and  these  did  no  better  than  simply  allowing  y  uncertainty. 

Fig.  lOd  shows  the  prediction  result  for  the  best-performing 
scenario  with  lowest  rms  wave-height  error,  highest  skill,  and  highest 
likelihood  ratio.  There  is  considerably  more  variability  in  the  wave- 
height  prediction  uncertainty,  compared  to  the  previous  examples. 
This  can  be  attributed  to  the  addition  of  the  wave  parameters  (y  and 
B ),  leading  to  an  increase  in  the  model  complexity  as  well  as  explicitly 
acknowledging  them  as  sources  of  uncertainty.  The  resulting  wave- 
height  prediction  uncertainty  was  higher  at  times  when  the  actual 
prediction  error  was  higher;  hence,  the  inclusion  of  model-parameter 
uncertainty  not  only  improved  the  prediction,  but  also  yielded  a  more 
accurate  error  estimate.  Specifically,  when  the  original  model  failed  to 
predict  the  largest  waves,  the  modified  model  either  improved  that 
prediction  (e.g.,  period  between  15  and  16  October)  or  updated  the 
uncertainty  such  that  it  included  the  observations  (e.g.,  the  period 
between  17  and  19  October).  Note  also  that  the  smallest  wave 
conditions  at  the  beginning  of  the  analysis  period  were  also  better 
predicted. 

The  wave-breaking  parameter  was  allowed  mild  uncertainty  such 
that  the  input  was  normally  distributed  with  a  mean  of  0.34  (the  prior) 
and  standard  deviation  of  0.05.  This  allowed  the  Bayesian  network  to 
attempt  to  predict  an  updated  parameter  value  based  on  the  offshore 
wave  height,  period,  direction,  and  water  depths.  Fig.  11  shows  the 
sensitivity  of  the  estimated  value  of  g  that  corresponded  to  the 
nearshore  wave-height  predictions  shown  in  Fig.  lOd.  Under  low 
wave  heights  (less  than  about  1.25  m),  the  most  probable  value  of  y  is 
about  0.3,  but  the  confidence  is  not  very  high  (probability  -50%)  and  a 
value  of  0.4  is  nearly  equally  likely.  However,  as  the  wave  height 
increased,  the  optimal  value  of  y  also  increased.  For  wave  heights 
exceeding  2  m,  the  optimal  value  of  y  is  0.4  or  higher.  And,  under  high 
waves,  the  parameter  uncertainty  decreased  (indicated  by  a  probability 
increase),  indicating  that  parameter  uncertainty  is  relatively  unimpor¬ 
tant  under  low  wave  heights  and  becomes  important  when  wave 
breaking  is  important.  This  dependence  on  wave  height  has  been 
demonstrated  by  others  (Apotsos  et  al„  2008;  Roelvink  and  Broker, 
1993;  Ruessink  et  al.,  2003).  Unlike  other  efforts  to  estimate  the  optimal 
parameter  values,  in  our  implementation,  the  prior  estimate  (y  =  0.34) 
continues  to  constrain  the  updated  estimate  because  we  do  not  allow 
large  uncertainties  in  this  parameter.  In  Part  II,  we  explore  the  problem 
of  inverse  modeling  and  the  impact  of  loosening  this  constraint. 


5.  Conclusions 

We  demonstrate  that  results  obtained  from  a  detailed  forward 
model  of  surf-zone  wave  evolution  are  reproduced  accurately  using  a 
Bayesian-network  modeling  approach.  We  tested  the  Bayesian 
network  using  both  simplified  scenarios  and  then  a  real-world 
prediction  example.  In  the  real-world  surf-zone  example,  prediction 
skill  was  about  80%  if  the  model  inputs  (offshore  wave  conditions, 
bathymetry,  and  free  parameters)  were  assumed  to  be  perfect.  Using 
both  simplified  scenarios  and  the  real-world  prediction  example,  we 
show  that  uncertainty  in  the  input  data  is  accurately  transferred  to 
uncertainty  in  the  predictions.  However,  if  uncertainties  in  the 
boundary  and  forcing  data  were  not  conveyed  to  the  Bayesian 
network,  overly  optimistic  prediction  uncertainties  were  computed. 
More  skillful  prediction  uncertainties  were  obtained  by  including 
measured  alongshore  bathymetric  variability  as  a  source  of  input 
uncertainty,  but  prediction  skill  under  this  scheme  improved  only 
marginally.  Alternatively,  if  model-parameter  uncertainty  is  included, 
the  prediction  skill  increased  substantially  (to  about  90%),  mean  and 
rms  errors  were  reduced,  and  the  predicted  uncertainties  were  more 
consistent  with  the  observations.  The  Bayesian  network  has  several 
advantages.  It  significantly  reduces  the  dimensionality  of  the  problem, 
compared  to  a  detailed  model;  uncertainty  estimates  are  made  for  all 
predictions,  and  it  is  possible  to  estimate  model  parameters 
simultaneously  with  making  the  wave-height  prediction. 
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